﻿using MathNet.Numerics.Interpolation;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using UtilZ.Dotnet.Ex.Base;
using UtilZ.Dotnet.Ex.Model;

namespace UtilZ.Dotnet.MathNetEx.Interpolation
{
    /// <summary>
    /// 线性插值扩展类,Piece-wise Linear Interpolation.
    /// </summary>
    /// <remarks>Supports both differentiation and integration.</remarks>
    unsafe public class LinearSplineEx : IInterpolation
    {
        private readonly int _length;
        private readonly double* _x;
        private readonly double* _c0;
        private double[] _c1;
        private readonly Lazy<double[]> _indefiniteIntegral;


        /// <summary>
        /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
        /// </summary>
        bool IInterpolation.SupportsDifferentiation => true;

        /// <summary>
        /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
        /// </summary>
        bool IInterpolation.SupportsIntegration => true;


        /// <summary>
        /// 构造函数初始化
        /// </summary>
        /// <param name="x">Sample points (N+1), sorted ascending</param>
        /// <param name="c0">Sample values (N or N+1) at the corresponding points; intercept, zero order coefficients</param>
        /// <param name="c1">Slopes (N) at the sample points (first order coefficients): N</param>
        private LinearSplineEx(SpanZ<double> x, SpanZ<double> c0, double[] c1)
        {
            this._length = x.Length;
            this._x = x.DataPtr;
            this._c0 = c0.DataPtr;
            this._c1 = c1;
            this._indefiniteIntegral = new Lazy<double[]>(ComputeIndefiniteIntegral);
        }


        /// <summary>
        /// Definite integral between points a and b.
        /// </summary>
        /// <param name="a">Left bound of the integration interval [a,b].</param>
        /// <param name="b">Right bound of the integration interval [a,b].</param>
        public double Integrate(double a, double b)
        {
            return this.Integrate(b) - this.Integrate(a);
        }

        /// <summary>
        /// Interpolate at point t.
        /// </summary>
        /// <param name="t">Point t to interpolate at.</param>
        /// <returns>Interpolated value x(t).</returns>
        public double Interpolate(double t)
        {
            int k = this.LeftSegmentIndex(t);
            return this._c0[k] + (t - this._x[k]) * this._c1[k];
        }




        /// <summary>
        /// Differentiate at point t.
        /// </summary>
        /// <param name="t">Point t to interpolate at.</param>
        /// <returns>Interpolated first derivative at point t.</returns>
        public double Differentiate(double t)
        {
            int k = this.LeftSegmentIndex(t);
            return this._c1[k];
        }

        /// <summary>
        /// Differentiate twice at point t.
        /// </summary>
        /// <param name="t">Point t to interpolate at.</param>
        /// <returns>Interpolated second derivative at point t.</returns>
        public double Differentiate2(double t)
        {
            return 0d;
        }


        /// <summary>
        /// Indefinite integral at point t.
        /// </summary>
        /// <param name="t">Point t to integrate at.</param>
        public double Integrate(double t)
        {
            int k = LeftSegmentIndex(t);
            var x = t - _x[k];
            return _indefiniteIntegral.Value[k] + x * (_c0[k] + x * _c1[k] / 2);
        }

        private double[] ComputeIndefiniteIntegral()
        {
            var integral = new double[this._c1.Length];
            for (int i = 0; i < integral.Length - 1; i++)
            {
                double w = this._x[i + 1] - this._x[i];
                integral[i + 1] = integral[i] + w * (this._c0[i] + w * this._c1[i] / 2);
            }

            return integral;
        }


        /// <summary>
        /// Find the index of the greatest sample point smaller than t,
        /// or the left index of the closest segment for extrapolation.
        /// </summary>
        private int LeftSegmentIndex(double value)
        {
            int index = SpanEx.BinarySearch(this._x, this._length, value);
            if (index < 0)
            {
                index = ~index - 1;
            }

            return Math.Min(Math.Max(index, 0), this._length - 2);
        }







        /// <summary>
        /// Create a linear spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x.
        /// </summary>
        public static LinearSplineEx InterpolateSorted(SpanZ<double> x, SpanZ<double> y)
        {
            if (x.Length != y.Length)
            {
                throw new ArgumentException("All vectors must have the same dimensionality.");
            }

            if (x.Length < 2)
            {
                throw new ArgumentException("The given array is too small. It must be at least 2 long.", nameof(x));
            }

            var c1 = new double[x.Length - 1];
            for (int i = 0; i < c1.Length; i++)
            {
                c1[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
            }

            return new LinearSplineEx(x, y, c1);
        }



    }
}
